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We investigate the mutual influence of impurities in two-dimensional d-wave superconductors involving self- 
consistent solutions of the Bogoliubov-de Gennes equations. The local order parameter suppression, the local 
density of states (LDOS) as well as the interference of impurity-induced structures are analyzed. We employ an 
impurity position averaging scheme for the DOS that does not neglect these interference effects, as the commonly 
used T-matrix approaches do. 
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In recent years scanning tunneling microscopy 
(STM) has become an excellent high resolution 
probe of the electronic structure near defects. In 
particular it has been used to investigate the ef- 
fects of individual Zn and Ni atoms in the d-wave 
superconductor BiaSraCaCuaOg+A- (BSCCO) |). 
These experiments show the appearance of a 
quasiparticle bound state with a fourfold symme- 
try around each impurity, reflecting directly the 
anisotropic d^2_y2 pairing symmetry [Q. 

The existence of impurity-induced bound states 
was predicted early on by Balatsky et al. using 
the T-matrix approximation |§,^. Often, local 
variations in the local order parameter (LOP) dy- 
namically generated by the impurity potential are 
neglected. Since the LOP varies only on a coher- 
ence length scale, which is very short (^ « 20 
A) in high-Tc superconductors, the overall ef- 
fect of this variation was argued to be negligi- 
ble. Recent calculations, however, indicate other- 
wise MH showing that the inclusion of the order 
parameter relaxation allows to achieve a better 
agreement with microwave conductivity data. 

In pure (1^2 _y2 superconductors the averaged 
density of states (DOS) p{E) vanishes linearly for 
\E\ — ^ 0. With impurities the self consistent T- 
matrix approximation m leads to a finite DOS 
at -E = 0, non perturbative approaches have var- 
iously predicted that p{E) vanishes Isj or that 
it diverges 0. These issues are still controver- 



sial and some efforts are needed to settle the low 
energy DOS issue in the presence of a finite con- 
centration of impurities. 

In this paper we include the local gap sup- 
pression of the order parameter and examine the 
LDOS for two nearby non-magnetic impurities 
and for a finite impurity concentration to ex- 
plore how impurity induced structures can inter- 
fere with each other, and to what extent this in- 
fiuences the impurity position averaged DOS. 

We start from a pairing Hamiltonian for elec- 
trons hopping on a square lattice with nearest 
neighbor hopping matrix element t, and a super- 
conducting bond order parameter A^^: 

{i,j),o- {1},<T 

The non-magnetic impurity potential is summed 
over impurity positions {{I}) and Ui is the local 
potential at the impurity site I . Energies will be 
measured in units of t and lengths in units of the 
lattice constant a. In what follows we set fj, = 
corresponding to half-filling since the variation of 
the filling does not change our conclusions. Note 
however that, even for ^ = 0, 7i is not particle- 
hole (PH) symmetric due to the presence of [/;. 
We first consider a single impurity located at 



(1) 



site i, and we choose C/; ~ Uq. We perform the 
unitary transformation [[lO[ 



Ci,( 



J2(ln,aU^,n + {-lVll. 
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and by imposing the condition that the transfor- 
mation (0) diagonaUzes the Hamihonian (nh we 
obtain the Bogoliubov-de Gennes (BdG) equa- 
tions [IToll 
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Here Aij is the self-consistent mean field solution 



for the order parameter 



A,j = V, 
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and Vi,j is the attractive pairing interaction. 

We solve numerically the self-consistent BdG 
equations (0) and (0) on lattices of 20x20 sites 
and with open boundary conditions. For the pair- 
ing interaction we choose 
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where Sij = 1 when i, j ^ i and i, j nearest neigh- 
bors. This form is motivated by the assumption 
that the pairing interaction is reduced around the 
impurity, and we therefore consider parameters 
with 1^1 1 < 1^0 1- Note that for Fi = Fq < this 
form of the pairing leads to a d^^-y^ structure of 
the order parameter. The local d-wave compo- 
nent of Ai ,- is obtained from: 



Ai — - iAi,i+x 
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In the clean case A^ is uniform in space and has 
a value of Af « 0.48 for Vb = -1- Modifications 
due to impurities or due to the border extend in 
space on a length scale of the order of the co- 
herence length ^. For our set of parameters we 
estimate ^ « 4 using ^ = vp/Ao, where vp is the 
Fermi velocity and Aq is the maximum of LOP. 
The value of Vq and consequently the value of Af 
are chosen to render ^ small enough on a 20x20 
system so that a sizeable portion of the lattice can 
still be considered as representative of the bulk of 
the system. Fig. |l| shows A^ near the impurity 
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Figure 1. Normalized d-wave component of the 
order parameter along the (1,0) direction for dif- 
ferent impurity potential strengths Uq. The pair- 
ing interaction is set to Vq = —1. Panels a,b,c 
correspond to different values of Vi to account for 
the suppression of the pairing interaction near the 
impurity. The position is measured with respect 
to the impurity site. 



site along the (1,0) direction for various values of 
Vi and Uq- A common feature is a strong sup- 
pression of the d-wave component of the order 
parameter around the impurity site on a range 
of roughly 3 lattice constants, in agreement with 
the estimated coherence length ^. The LOP along 
the (1,1) direction is modified on the same length 
scale. Note that in the case of Vi = we have 
Af = since Ai,;±:r = Ai;,±y = 0. 
The LDOS as obtained from 

p.iE) = ^ [K\'S{E - En) + K\'S{E + En)] (7) 

n 

near the impurity is plotted in Fig. y| as the con- 
tinuous line. In Eq. (|^) the sum is over the eigen- 
states n of the Hamiltonian (Q) and En denote 
the corresponding eigenenergies. In the numeri- 
cal calculation the delta functions present in Eq. 



([7|) are approximated by Lorentzians with a width 
r. Typically we use T/t = 0.1. On the impurity 
site the LDOS is strongly reduced while on the 
nearest neighbor sites we find two peaks around 
zero energy in the LDOS, associated with two 
resonant quasiparticle bound states with energies 
±ujR, (see Fig. ^ where the DOS is calculated 
with r = 0.001). The two peaks are not resolved 
in Fig. H due to the large broadening F = 0.1, 
and only one peak is observed around zero en- 
ergy. This peak is, also if not evident from the fig- 
ure, asymmetric revealing the underlying double 
peak structure at :\iujr. Measuring the distance 
variation of the LDOS around the impurity at 
the resonance energy reveals also the well known 
fourfold symmetric structure characteristic for a 
dx2_y2 superconductor ||]. 

In the next step we add a second impurity 
and explore the mutual influence of the two im- 
purities. The resulting LDOS at a selected site 
between the two impurities is shown in Fig. y 
(dashed line), for different impurity distances 
along the (1,0) direction. With one impurity we 
observe a resonance peak at very low energy oju 
corresponding to the lowest positive eigenvalues 
of the system that has an energy of Eq — loh « 
0.04. When a second impurity is added at a dis- 
tance d = 2 from the first, the peak is strongly 
enhanced by a factor of roughly 2.5. Its asym- 
metry is again reflecting the underlying resonant 
impurity states at positive and negative energies 
(see Fig. ||) . When the second impurity is moved 
one site further away, the peak disappears leav- 
ing a dip around zero energy. The peak is again 
recovered when the distance is increased further 
(d = 4). This alternating variation is naturally 
explained as an interference phenomenon due to 
the overlap of the bound state wave functions 
near each impurity. In order to demonstrate in 
detail what is happening at low energies when 
one or more impurities are added to the system 
we plot in Fig. | the DOS for F = 10"^ In Fig. 
g|.a we show the clean case as a reference; the 
spikes at \E\ k, 0.1 are the lowest energy eigen- 
values for a 20x20 system with open boundary 
conditions; note the symmetry of the LDOS in 
the clean case: for Uq — Q the Hamiltonian is 
particle-hole symmetric for /i = 0. When we add 
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Figure 2. Local density of states in the presence of 
two impurities measured near one impurity site. 
The continuous lines are the calculated densities 
of states with only one impurity, the dashed ones 
when the second impurity is added. The LDOS 
is always measured immediately on the adjacent 
site of the first impurity. The parameters used 
are: F = 0.1, Vq = -1, Vx = 0, t/o = 8. 



an impurity at the center of the system we observe 
the appearance of two additional eigenvalues with 
asymmetric weights at energies \E\ « 0.04. The 
asymmetry results from the broken particle-hole 
symmetry due to the presence of the impurity. 
When we add a second impurity near the first 
one (see Fig. ^c) the number of eigenvalues in 
the region \E\ < 0.1 is doubled. Thus, by adding 
impurities one by one an impurity band of finite 
width is formed due to the overlaps of individual 
impurity induced quasiparticle bound states [O . 
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Figure 3. Plot of the DOS with F = 0.001 in 
the region \E\ < 0.1. a) Clean case result; b) 
with one impurity; c) with two impurities located 
three sites away one from each other in the (1,0) 
direction, for Uq = 8, Vq = —1, Vi = 0. 



When a self-consistent T-matrix approach is 
used to calculate the impurity position averaged 
DOS the interaction effects between impurities 
are neglected due to the selection of non-crossing 
diagrams for impurity scattering in the standard 
T-matrix approximation. Therefore, the mu- 
tual influence of impurities as observed above is 
missed in the conventional T-matrix approach. 
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An obvious question to ask is to what extent 
the interference effects and the impurity band 



Figure 4. Averaged DOS for an impurity concen- 
tration of rii = 20%. The inset is a zoom of the 
low energy part of the calculated DOS. 



formation may survive an appropriate impurity 
position averaging scheme. In order to answer 
this question we introduce in the system a fixed 
number of impurities in random positions and 
calculate the LDOS in the center of the lattice 
P^^\E). We keep the site for the measurement 
of the LDOS fixed for all impurity configurations 
with the same impurity concentration. The av- 
eraged DOS is calculated according to: where 
m enumerates an impurity configuration and Nc 
is the total number of impurity configurations. 
Typically we choose Nc ^ 10^, convergence is 
usually achieved already for Nc > 50. In Fig. 
H Pav{E) is shown for a large impurity concen- 
tration rii — 20% to highlight a new low energy 
structure in the averaged DOS that is not ap- 
pearing in the usual T-matrix calculation, that 
leads to a flat and finite DOS around zero en- 
ergy. The new feature that we find consists of an 
asymmetric structure centered around zero en- 
ergy. This previously unresolved feature in the 
DOS survives also when n^ decreases to a few 
percent although it becomes less pronounced; its 
origin is clearly related to the overlap of the im- 
purities induced bound states and the accompa- 
nying impurity band formation. The above dis- 
cussed evolution of the impurity band naturally 
suggests furthermore that the DOS at zero en- 



ergy will remain zero even for a finite impurity 
concentration. 

In conclusion we have outlined an extension of 
the BdG formalism to two dimensional supercon- 
ductors with a finite concentration of impurities. 
We have shown how impurities interact between 
each other giving rise to interference phenomena 
and leading eventually to an impurity band for- 
mation due to the finite overlap of quasiparticle 
bound states. These phenomena survive for an 
appropriately chosen impurity position averaging 
scheme leading to a double peak structure in the 
averaged DOS at low energies which may be re- 
solvable in STM measurements. 
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